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Surface irregularities such as protuberances inside a hypersonic boundary layer may lead to premature 
transition on the vehicle surface. Early transition in turn causes large localized surface heating that could 
damage the thermal protection system. Experimental measurements as well as numerical computations 
aimed at building a knowledge base for transition Reynolds numbers with respect to different protuberance 
sizes and locations have been actively pursued in recent years. This paper computationally investigates the 
unsteady wake development behind large isolated cylindrical roughness elements and the scaled wind-tunnel 
model of the trip used in a recent flight measurement during the reentry of space shuttle Discovery. An 
unstructured mesh, compressible flow solver based on the space-time conservation element, solution element 
(CESE) method is used to perform time-accurate Navier-Stokes calculations for the flow past a roughness 
element under several wind-tunnel conditions. For a cylindrical roughness element with a height to the 
boundary-layer thickness ratio from 0.8 to 2.5, the wake flow is characterized by a mushroom-shaped 
centerline streak and horse-shoe vortices. While time-accurate solutions converged to a steady-state for a 
ratio of 0.8, strong flow unsteadiness is present for a ratio of 1.3 and 2.5. Instability waves marked by distinct 
disturbance frequencies were found in the latter two cases. Both the centerline streak and the horse-shoe 
vortices become unstable downstream. The oscillatory vortices eventually reach an early breakdown stage 
for the largest roughness element. Spectral analyses in conjunction with the computed root mean square 
variations suggest that the source of the unsteadiness and instability waves in the wake region may be traced 
back to possible absolute instability in the front-side separation region. 

I. Introduction 

H ypersonic Boundary layer transition is one of the active research topics in recent years due to its importance in 
many aerodynamic and aerothermodynamic applications. Laminar-turbulent transition onset and the 
accompanying rapid rise in surface heating can directly impact the performance and safety of hypersonic vehicles. 
From the design standpoint, accurately predicting transition may result in a much more efficient thermal protection 
system (TPS). In addition to natural transition in response to environmental perturbations, large surface 
irregularities such as protuberances or cavities are known to lead to premature transition on the vehicle surface. 
Large size protuberances thus can potentially cause severe damage due to high surface heating rate. The repair of 
protruding gap fillers by an astronaut to prevent TPS failure during the Shuttle mission STS- 114 is one good 
example for such consideration. Due to the lack of a robust and efficient computational tool for flowfield prediction, 
empirical relations based on wind tunnel and flight data are used to estimate transition in the presence of surface 
roughness. In order to enhancing current understanding of the flow physics and physics based prediction tools for 
future hypersonic vehicle design, a combination of experimental measurements and higher fidelity numerical 
computations of viscous flow over large surface irregularities constitutes one of the most important goals in 
hypersonic research. 

Most of the computational investigations for the hypersonic roughness element problem use structured mesh Navier- 
Stokes equations solvers. To handle the protruding geometry, researchers have used the immersed boundary 
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formulation in conjunction with existing time-accurate compressible Navier-Stokes solvers to calculate the wake 
flowfield [1-4]. The roughness height being simulated is relatively small in these computations due to the limitation 
in handling complex shock patterns of the underlying high-order schemes embedded in these Navier-Stokes solvers. 
In principle, an unstructured mesh is preferred because it offers more flexibility in grid generation and can 
potentially save computational time by avoiding unnecessary grid clustering due to topological constraints common 
in structured meshes. To handle hypersonic viscous flow over a protruding roughness element in the boundary 
layer, several known numerical issues must be addressed. Most existing unstructured mesh CFD codes are based on 
well-known upwind schemes that utilize a ID approximate Riemann solver at cell interfaces. Unique flow 
dependent variables determined by the Riemann solver are required at these cell interfaces to integrate the governing 
conservation laws over the discretized domain. When applied to multi -dimensional flows, such one-dimensional 
approximations give rise to undesired numerical problems, especially for triangular or tetrahedral unstructured 
meshes [5]. For hypersonic applications, the carbuncle phenomenon takes place near the blunt nose region [6] and 
accurate computation of viscous flows with tetrahedral meshes remains a challenge [7]. Remedies for these 
numerical issues have been proposed and investigated but in general, lack robustness. 

When a large surface roughness is submerged in a hypersonic boundary layer, complex flow physics such as shock- 
shock, shock-boundary layer interactions, flow separation, unsteady shear layers and wakes may develop around and 
downstream of the roughness site. In an attempt to resolve such complex flow physics, a time-accurate, 
unstructured mesh Navier-Stokes solver, ez4d , [8] based on the space-time conservation element, solution element 
(CESE) method [9] was used in our previous study to investigate hypersonic flows over a rectangular prism and 
cylindrical roughness element [10]. Several configurations from NASA Langley wind tunnel experiments [11] were 
computed. It was shown that at a boundary layer edge Mach number around 4, the wake region becomes unsteady 
for a rectangular fence with a height of about 1.3 times the boundary-layer thickness. Unlike its low speed 
counterpart, strong vortex shedding was not observed for both rectangular fence and cylindrical roughness elements 
at a Mach number of 4 or 6.5. Relatively coarse hexahedral meshes were used for the cylindrical roughness element 
during this study. 

The main objective of this paper is to continue the 3D Navier-Stokes computations for the cylindrical and the 
Detailed Test Objective (DTO) roughness trips investigated in NASA Langley’s wind tunnel experiments [12]. The 
same cylindrical roughness element with a height of 2 mm and a diameter of 4mm as investigated in [10] is studied 
with a more refined tetrahedral mesh and in addition, for a higher Reynolds number. Three wind tunnel conditions 
for the cylindrical roughness element are investigated in detail to help build the knowledge base concerning the 
effects of the roughness height and roughness Reynolds number. Flow visualization and unsteady data processing is 
used to explore detailed flow physics and possible instability wave development in the wake region. The time- 
accurate solutions are also analyzed spectrally to understand the nature of the instability waves in the roughness 
wake. In the following section, a brief description of the numerical method used is given, followed by results and a 
brief summary. 


II. Numerical Method 


Three-dimensional compressible Navier-Stokes equations in vector form can be written as 


8Q 8F 8G 8H „ 8F V 8G V 8H V 

8t 8x 8y 8z 8x 8y 8z 


( 1 ) 


where the dependent variable vector is defined as Q = (p, pu, pv, pw, e) and v, y, z, and t represent spatial 
coordinates and time, respectively. Flow variables p, u, v, w, and e represent density, three velocity components, and 

total energy (e = p l(y — X) + p(u 2 + V 2 + W 2 )/ 2 , where p denotes the pressure), respectively. Definitions of 
the inviscid flux vectors F, G, H and viscous flux vectors F v , G v , H v can be found in standard text books and will not 
be included here. The source vector S contains all external forcing or other energy -related source terms. To close the 

system, the perfect gas relation, p = pRT , with T representing the temperature, is used in conjunction with eq. (1). 
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The space-time CESE method is formulated in the strong form of the flux equations. The space-time flux h is 
defined as: 


h=(F-F v ,G-G v ,H-H v ,Q ) (2) 

By using Gauss’ divergence theorem in the space-time domain, eq. (1) is rewritten in the following integral form 

| h-ds+^SdV =0 (3) 

v 

where the space-time flux vector is integrated over the surface T of an arbitrary space-time volume V. The surface 
normal vector is defined by s = n dA , where dA is the area increment on T, and n is the outward unit normal 

vector. Equation (3) is quite general, and in fact, any conservation law with or without a source vector can be cast in 
this form. Thus, the numerical algorithm devised for eq. (3) can be easily extended to other physical problems. 

The space time CESE method is formulated by enforcing flux conservation over a discretized space-time domain. 
This can be illustrated with the triangular elements in a two-dimensional space shown in Fig. 1 . Numerical solution 
of the dependent variables within the triangular element BDF is assumed to satisfy the following first order 
equation: 


Q(x,y,t) = Q 0 +Q,(t-t 0 ) + Q x (x-x 0 ) + Q y (y-y 0 ) (4) 

where Q 0 is the solution vector at the solution point O. The vectors Q t , Q x , Q y , and Q z are analogous to the 
derivatives of the dependent variables. In contrast to the conventional unstructured mesh method, flux integration is 
carried out not on the triangular element directly. Instead, three conservation elements (CE) defined by the 
quadrilaterals GBCD , GDEF , and GFAB around three neighbors are used for the integration of the conservation 
laws. In the space-time domain, these three quadrilateral CEs extend to three quadrilateral prisms. The solution 
point O is defined as the geometric center of three conservation elements, ABCDEF. The discretized flux equation 
for all three prisms can be simplified to 

(A + A + A)Qo + + h + ^v (5) 

where A h A 2 , and A 3 are the surface areas of GBAF , GFED , and GDCB, respectively. Discretized integrals I b / 2 , 
and / 3 are flux integrals for outer-side and bottom faces for three surrounding elements. Source terms in eq. (3) can 
be integrated to form S v . Semi-implicit treatment of this source term is usually done to improve solution accuracy 
and avoid stiffness issues. Equation (5) is an algebraic equation and can be solved easily without any matrix 
inversion. The use of eq. (4) results in a formally second-order accurate scheme. Discretized equations including all 
geometric details for the Euler equations can also be found in [13]. 

In the above discretized equation, flux integrations are carried out on all faces of the conservation elements shown in 
Fig. 1(b). Flux vectors are evaluated at faces AB, BC, CD, DE, EE , and EA in the integrals I b / 2 , and I 3 These faces 
reside over the smooth regions of all three neighboring elements. As a result, flux vectors are uniquely defined in 
the discretized equations, thus no exact or approximate Riemann solver is needed. In light of all the numerical 
issues associated with multi-dimensional upwind methods, this distinct feature of the CESE method offers a 
numerical framework that is superior to other methods in terms of the flux conservation properties in the discretized 
domain. 

The derivatives Q t , Q x , and Q y , in Eq. (4) are yet to be determined. In the non-dissipative a-scheme, these 
derivatives are computed by solving two remaining flux equations (in addition to the summation of all three 
conservation elements, Eq.(5)) in conjunction with the following equation: 

Q t +AQ x + BQ y =S 
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where A and B are the Jacobian matrices of the flux vectors F and G, respectively. Thus, there are a total of three 
equations for three unknowns. Unfortunately, the a-scheme is only neutrally stable for the nonlinear Euler and 
Navier-Stokes equations and cannot be used in general. Nevertheless, the formulation of such a dissipation-free 
scheme serves as a baseline to gauge how much numerical dissipation is being added in the dissipative scheme used 
in the Euler or Navier-Stokes equations. In practice, derivatives are evaluated by using solution variations around 
three neighbors. In the original CESE method, three neighboring solution points are used to obtain three sets of Q x , 
and Q y . The unique values are determined by applying a weighted average over three sets of derivatives. The 
weighted average helps to eliminate oscillations around flow discontinuities. While this approach works well for 
relatively uniform meshes, numerical problems arise for highly non-uniform meshes where the triangles used for 
derivative evaluation are inconsistent with the CE faces. The edge-based derivative approach suggested in Ref. [8] 
modifies the derivative evaluation procedure to remedy the above shortcomings. In this approach, the three vertices 
of the triangular elements with dependent variables evaluated by using three neighboring solution elements are used 
to obtain three sets of neighboring derivatives. The same weighted average procedure of the original CESE 
approach is then used to uniquely determine the derivatives. There are several weighted average procedures 
suggested by Chang [14]. In this paper, the second weighting function described in Ref. [14] is found to be the most 
effective and robust. The weighted average is analogous to flux limiters used in conventional upwind schemes. 
However, it should be noted that in the CESE method, one universal weighted average works for all problems. No 
tuning for shocks or stagnation regions is necessary. 




Figure 1. Conservation elements for a 2D triangular mesh: (a) top view (b) 
space time CE volume 


The above derivation is valid for 2D triangular elements but it can be extended easily for 3D equations. There are 
questions often raised concerning the fact that a picture of 3D conservation elements similar to that shown in Fig. 1 
for 2D cases cannot be drawn. In practice, the CE manifolds need not to be drawn for 3D solutions. As long as the 
physical time is uniquely defined on each CE face in the flux integration procedure, it is irrelevant how these CE 
manifolds look in a four-dimensional space. Both hexahedral and tetrahedral meshes are used in the present study 
for hypersonic viscous flow computations. 

The other issue often raised concerns the fact that for the CESE method, numerical dissipation increases 
significantly when the CFL number is decreased. For this reason, steady state calculations are done by using local 
time-stepping, i.e., constant CFL number for all elements. Normally, a CFL number as close to one as possible is 
used to ensure accuracy. For time-accurate calculations, a constant time step implies that the CFL number varies 
substantially over the domain if a highly non-uniform mesh is used. The Courant number insensitive (CNI) method 
[13] should be used to avoid excessive numerical dissipation in the coarse mesh region. For the edge -based 
derivative approach, the vertices used for derivative calculations are pushed away from the solution point to increase 
numerical dissipation and pushed the opposite direction to reduce the dissipation. A parameter a is used to control 
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numerical dissipation. A a value of 1 means zero dissipation. In a CNI scheme, the a value is very close to 1 in the 
coarse grid region and increases to a larger value (around 1.3-1. 5 for subsonic flows, and 5 or higher for supersonic 
flows with shocks) to maintain numerical stability in the highly stretched mesh inside the boundary layer. 

III. Results and Discussion 

Hypersonic flows over an isolated cylindrical roughness and the DTO trip are the main focus of the present paper. 
Three flow conditions taken from NASA Langley Research Center’s Mach 10 wind tunnel experiments by Danehy 
et al. [11] were computed for both trips. The schematic of the experimental setup is shown in Fig. 2 for a cylindrical 
roughness element. The free-stream Mach number is 9.65 and the free-stream temperature is 53.3 K. The 
conventional Sutherland law is used throughout the computations presented herein. Since a different viscosity law 
was used to compute the unit Reynolds number reported for the experiment [11], a modified value is used to mimic 
one of the experimental conditions (specifically, condition II defined below) near the roughness element. The free- 
stream unit Reynolds number for conditions I and II is 3.513 x 10 6 per meter, which is calculated based on an 
estimate of the post-shock Reynolds number from the experiment [11]. The free-stream unit Reynolds number for 
condition III, 7.4 x 10 6 /m, corresponds to the largest Reynolds number in the wind tunnel experiment with a 
stagnation pressure of 1350 psi. The wall temperature is 308 K. In all the experiments, the roughness element is 
mounted on a flat plate with a free-stream flow incident angle of 10° or 20° for Conditions I and II-III, respectively. 
The corresponding post-shock Mach numbers are 6.52 and 4.16 for an angle of attack of 10° or 20°, respectively. 
The cylindrical roughness has a 4 mm diameter and 2 mm height. The DTO trip geometry to be described later has a 
height around 1 mm. The ratio of roughness height to boundary layer thickness is estimated to be about 0.82, 1.3, 
and 2.5 under Conditions I, II, and III, respectively. Such roughness height is relevant for many hypersonic 
applications. The length Reynolds number is about 3.84 x 10 5 , 2.6 x 10 5 , and 5.4 x 10 5 at the roughness center. The 
roughness elements are either in the subcritical or slightly larger than the critical Reynolds number at which the 
laminar flow becomes unstable to first- or second-mode instability waves. Table 1 summarizes relevant parameters 
for the cylindrical roughness under all three conditions in which M e indicates the post shock Mach number, k/8 the 
roughness height to boundary layer thickness ratio, Re x the Reynolds number at the roughness location, Re h the 
roughness height Reynolds number based on the post-shock conditions, and Re k the roughness Reynolds number 
evaluated by using the local flow conditions at the roughness height. For all computations shown in this paper, only 
half of the domain is computed. Symmetry conditions are applied at both span wise boundaries. In all the 
computations reported in this paper, the inflow boundary is located 75.4 mm upstream of the roughness center (or 
approximately 72.4 mm from the leading edge for experiments with the cylindrical roughness). At the inflow, a 
converged solution obtained by a 2D Navier-Stokes calculation is imposed. Constant temperature viscous boundary 
conditions and non-reflecting conditions are employed at the wall and the outflow boundaries, respectively. 

Time-accurate Simulations of an Isolated Cylindrical Roughness Element 

In Ref. [10], hypersonic flows over a cylindrical roughness element were investigated for conditions I and II using 
hexahedral meshes. For better solution accuracy, a tetrahedral mesh is preferred. Grid refinement studies were 
perfomed using several tetrahedral meshes with a mesh size ranging from 8 to 33 million elements. The finest 
tetrahedral mesh shown in Fig. 3 with about 33.1 million tetrahedrons was found to give the best resolution for both 
the horse-shoe vortices and the wake streak for all three flow conditions although the second finest mesh with about 
16.5 tetrahedrons gave very close solutions for condition I. In general, it is difficult to judge the current solution is 
already grid-converged without an extremely expensive grid-convergence study. Compared to the hexahedral mesh 
used in [10], the current tetrahedral mesh has better grid resolution in both stream wise and span wise directions. The 
spanwise grid resolution appears to be relatively more important to resolve both the steady and unsteady features of 
the flowfield. Figure 4 compares the streamwise evolution of velocity contours on the cross planes at four locations 
downstream of the roughness under condition I obtained by using the 5.6 million hexahedral elements and the 
current 33.1 million element tetrahedral mesh. Substantial improvement in the flow features is evident. 
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Figure 2. Schematic of the isolated cylindrical roughness element investigated in the LaRC wind tunnel 

experiments 


Table 1 Three free-stream conditions investigated for the isolated cylindrical roughness and the 
corresponding Reynolds number based on the height and locations 


Flow Condition 

M e 

k/S 

Re x 

Re h 

Re k 

I 

6.52 

0.82 

3.84 x 10 5 

10,600 

5,940 

II 

4.16 

1.3 

2.6 x 10 5 

6,800 

6,130 

III 

4.16 

2.5 

5.4 x 10 5 

14,300 

13,300 



Figure 3. Tetrahedral mesh with 33.1 million elements used for cylindrical roughness element 
computations: (a) 3D domain (b) symmetry plane and wall mesh (c) close up of the roughness element 
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For all time-accurate computations shown here, local time stepping was first employed for 30,000 iterations to have 
a roughly settled solution. Time-accurate computations were then carried out until the solution either converged to a 
steady state or reached a quasi-steady state for unsteady flows. Figures 5-7 show time traces of stream wise velocity 
and temperature at five probe points defined in Table 2 obtained from the time-accurate computations under three 
flow conditions. Condition I apparently reaches a steady state while both conditions II and III remain unsteady 
throughout the time-accurate calculations. However, solutions at all five probes have reached a large time 
asymptotic state for the latter two cases since the oscillatory time traces are repeated multiple times. A dominant 
frequency can be discerned in condition II. In contrast, multiple -frequency waves are present for condition III. 
These unsteady time traces suggest possible flow instability in the wake region which will be discussed in more 
detail later. 



Figure 4. Comparison of streamwise velocity contours at four downstream locations for steady-state solutions 
under condition I using two different meshes: (a) 5.6 million hexahedral mesh (b) 33.1 million tetrahedral 
mesh 

Table 2. Coordinates and locations of the probe points for the cylindrical roughness 


Probe 

X 

y 

z 

Location 

PI 

-2.13 

1.00 

1.90 

Front side separation region 

P2 

3.16 

0.92 

0.16 

Right behind the roughness 

P3 

11.18 

1.68 

0.16 

Near wake 

P4 

30.83 

1.57 

3.90 

Middle wake 

P5 

44.35 

2.97 

0.16 

Far wake 


Details of the solution can be better visualized by the numerical Schlieren pictures at the symmetry, the exit, and the 
y = 1.5 mm planes shown in Figs. 8-9. The leading edge shock, front-side separation region induced shock and the 
wake region recompression shocks are clearly visible. The signature of the horse-shoe vortex is also visible at the 
exit plane. The bright white curve at about the same height of the roughness in the wake region in Fig. 8 marks the 
boundary of the streak right behind the cylinder. For conditions II and III, these streaks are unsteady and oscillatory 
along the streamwise direction. For condition III, additional small structures develop around the oscillatory streak. 
The development of the horse-shoe vortices around the cylinder cross-section can be clearly visualized in Fig. 9. 
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For the steady vortices under condition I, these vortices remain straight throughout the computational domain. In 
contrast, horse-shoe vortices meander downstream for conditions II and III. It is interesting to note that the streak 
has a span wise symmetric, vertical-swinging oscillatory motion (a natural outcome of the symmetric boundary 
conditions at the centerline) while the horse-shoe vortices oscillate in a spanwise side-swinging fashion. While both 
of the larger roughness cases give rise to wake instability, the flow under condition III is apparently much more 
unstable because the oscillatory motion is much more pronounced and some additional structures are present near 
the exit (see Fig. 9(c)). 

The resulting surface limiting streamlines at a selected instant of time are plotted in Fig. 10 for all three conditions 
along with the vorticity magnitude contours at the exit plane. The limiting streamline patterns, similar to that 
observed in [10], show substantial upstream influences in the front-side separation region. For condition III, 
limiting streamlines appear to be wavier than the other two conditions. Solutions presented here are instantaneous 
snapshots, time-average behavior on the surface may not show any waviness. Again, both the symmetry plane 
streak and horse-shoe vortices are more pronounced at the exit plane in conditions II and III. Figure 1 1 depicts the 
normalized surface temperature gradient (normalized by the free -stream temperature and the length unit of 0.001 m) 
distributions. The images from Fig. 1 1 further confirm the existence of a streaky pattern around the cylinder with 
more streaks for relatively larger roughness elements. These streaks are related to vortex formation around the 
cylinder. More streaks are present for relatively more unstable trips. This behavior is in line with the surface 
thermography observations for the DTO trip by Danehy et al. [12]. 




(a) (b) 

Figure 5. Time histories at five probe points defined in Table 2 for cylindrical roughness element under 
Condition I: (a) streamwise velocity (b) temperature 

To visualize the overall vortex pattern around the roughness, the 3D isosurfaces of the vorticity magnitude are 
shown in Fig. 12. The vortex pattern in Fig. 12(a) shows a stationary, and fairly straight shape horse-shoe and wake 
streak vortex tubes for condition I. In contrast, the vortex tube formed by constant vorticity magnitude meanders 
downstream for both conditions II and III. For the most unstable case, condition III, the vortex evidently breaks 
down near the end of the computational domain. Smaller structures are generated when the breakdown takes place. 
More grid refinement studies are necessary to determine whether the vortex breakdown process was properly 
resolved in the computations. Figure 13 shows the slices of vorticity magnitude contours at a horizontal plane, y = 
0.5 mm. For comparison, the PLIF image from Ref. [11] for the flow at a similar height under condition III is also 
included in the figure. Similar large disturbances are evident both in computations (Fig. 13(c)) and experiments 
(Fig. 13(d)). It is interesting to note that more streaks upstream and slightly downstream the cylindrical roughness 
element are present as the flow becomes relatively more unstable. These streaks are also related to the streaky 
patterns in surface heat transfer rate shown in Fig. 11. The streamwise evolution of the vortex pattern can also be 
viewed from Figs. 14-16 where the vorticity magnitude contours are plotted at four streamwise locations with x = 5, 
25, 45, 65 mm. All coordinate values (x and z ) present herein are in mm. The stable case, condition I, has a stable 
and slowly rising horse-shoe vortex and streak. On the other hand, both the streak and the horse-shoe vortices rise 
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up substantially from the wall for conditions II and III. Signs of horse-shoe vortex breakdown manifested by 
weakened high shear layers on top of the horse-shoe vortices at v = 65 mm are also evident for condition III. 




Figure 6. Time histories at five probe points defined in Table 2 for cylindrical roughness element under 
Condition II: (a) streamwise velocity (b) temperature 




(a) (b) 

Figure 7. Time histories at five probe points defined in Table 2 for cylindrical roughness element under 
Condition III: (a) streamwise velocity (b) temperature 


Stability Analysis of the Flowfield around the Cylindrical Roughness Element 

Time traces at the probe locations and the vortex structures described above show strong evidence of instability 
waves for conditions II and III. Note that these instability waves are self-excited by numerical noise. Work is 
underway to investigate controlled disturbances input. In this paper, some preliminary stability analyses for the 
computed flow under condition II were performed to help understand the nature of the self-excited instability waves 
in the cylindrical roughness wake. During the time-accurate computations, the time-average of all dependent 
variables and their squares were computed and stored. These time-averaged solutions allow stability analyses of the 
mean flowfield. Flow disturbances can also be computed by subtracting the mean values from the instantaneous 
ones. The root mean square (rms) of the flow disturbances can be obtained by the computed time-averaged values 
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of square quantities such as u, T 2 , etc. Figure 17 shows the time-average stream wise velocity contours at v = 40 
mm. As shown in Refs. [4, 16], a well-formed streak and horse-shoe vortex pattern is expected to give rise to 
instability mechanisms that are similar to crossflow or Gortler vortices. Highly inflectional profiles in these vortical 
structures are susceptible to inviscid instability. Similar partial differential equation (PDE) based, 2D eigenvalue 
analysis has been performed for the meanflow shown in Fig. 17 to analyze possible instability modes of these wake 
structures. The methodology used for this analysis is described in Ref. [16]. For such meanflow, instability waves 
for the streak and horse-shoe vortex have different instability wave characteristics. Yet, it is non-trivial to separate 
these two flow structures in the 2D eigenvalue analysis. To isolate the instability of the streak, the span wise domain 
is truncated to z- 4 mm in the analyses. There are 32 span wise Fourier modes and 121 wall-normal points used in 
the calculations. Only symmetrical modes are analyzed because the time-accurate solutions shown above suggest 
the existence of a symmetrical instability for the streak. The stability results are shown in Fig. 18. Temporal 
analyses were performed and then transformed to spatial growth rates using the well-known Gaster transformation 
[17]. Common in flows with highly inflectional profiles, there exist multiple unstable modes with the most unstable 
modes around 90 and 150 kHz (modes 1 and 2 shown in Fig. 18(a)). The eigenfunction structures are concentrated 
in the high shear region on the top and sides of the mushroom-like streak structure. Isolating the horse-shoe vortex 
is more difficult. The domain was truncated at z— 12 mm with most points clustered around the horse-shoe vortex 
between z = 4 and z = 7 mm. Only anti -symmetric modes were investigated because the time-accurate simulations 
above suggest a spanwise side-swinging mode. The computed results are shown in Fig. 19. Due to the relatively 
low resolution, only one good-quality mode was tracked. The most unstable frequency is around 215 kHz. The 
eigenfunction structures for this and a lower frequency 41 kHz mode are also shown in the figure. Most 
perturbations are again in the high shear region. 

To facilitate comparison with the above stability results, time traces of stream wise velocity and temperature at a 
probe point P6 = (41.6 mm, 1.9 mm, 5.9 mm) are analyzed spectrally using the discrete Fourier transform. Figure 
20 shows the time traces and the resulting Fourier spectrum of temperature perturbations at both P6 and P5 (defined 
in Table 2). The Fourier mode amplitude shown in the figure has been normalized by the corresponding mean 
values. Point P6 is located very close to the x = 40 mm plane and P5 is further downstream but well within the 
wake streak. At both locations, there is a clear spectral peak at around 40 kHz. The second and third peaks are 
around 75 and 115 kHz, respectively. In fact, the 40 kHz mode is present and dominant throughout most of the 
wake region and right in front of the cylinder. The second and third peaks are quite close to the first and second 
harmonic of the primary 40 kHz mode. This self-excited disturbance frequency is lower than that predicted by the 
2D eigenvalue analysis. However, the eigenvalue analysis results were only valid at one specific stream wise 
location. To track the most unstable mode, eigenvalue analysis must be carried out at consecutive streamwise 
locations starting from slightly upstream of the roughness, which is quite computationally demanding and beyond 
the scope of this paper. Refs. [3-4] demonstrate the application of integrated amplification characteristics of wake 
instabilities for the purpose of analyzing transition due to periodic and isolated roughness elements, respectively. 
The computed rms streamwise velocity and temperature amplitude contours at the x = 40 mm plane are shown in 
Fig. 21. The self-excited disturbances are mainly concentrated in the horse-shoe vortex region, with a weaker 
perturbation sitting on top of the streak. The rms amplitude contours bear some similarity to those obtained from the 
eigenvalue analyses shown in Figs. 18-19. 

Results from the above secondary instability analysis at a representative location downstream of the cylindrical 
roughness above suggest that the most unstable instability mode in the wake region may have a frequency higher 
than the 40 kHz mode observed in the simulations. But then why is this relatively low frequency mode being 
excited by numerical noise? To answer this question, time traces at the probe point PI was further analyzed 
spectrally and plotted in Fig. 22. Point PI is located slightly upstream of the cylinder within the front-side 
separation region. There is a clear spectral peak at 40 kHz for both streamwise velocity and temperature at this 
location. Although not shown here, other probe points nearby also have an identical peak frequency. The 
disturbance structure at this dominant frequency can be visualized by the rms amplitude of the streamwise velocity 
plot at a nearby streamwise plane, x = -2 mm, as shown in Fig. 23. Large disturbances reside in the vicinity of high- 
shear regions roughly on top of the negative streamwise velocity separation bubble. The separation region on this 
streamwise plane results in a doubly inflectional profile (with streamwise velocity changing signs twice) around z = 
3-4 mm. Such profiles may (or may not) cause an absolute instability observed in separation regions at low-speed 
boundary layers. The existence of a distinct disturbance frequency around and after the cylinder and only one 
distinct spectral peak seems to lead to a likelihood of absolute instability. To trace the origin of this mode further 
upstream, the rms amplitudes of the temperature fluctuations at all grid points on the symmetry plane (z = 0) are 
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plotted in Fig. 24. In such a plot, values at different wall-normal distances are shown as discrete points. Tracking 
precise evolution of a particular point is difficult but the plot allows inspection of the variation along the streamwise 
direction. As can be seen, there are two distinct sharp peaks around x = -4 mm and x = -2 mm, respectively. After 
these spikes, temperature disturbances decay until the second bout of growth takes place around v = 10 mm. Note 
that the growth rate downstream of the wake appears to be much smaller than the two sharp peaks in the front-side 
separation region. The fact that this highly localized, rapid spatial growth occurs within the reverse flow region and 
corresponds to an almost constant frequency further suggests the existence of absolute instability. If this argument 
holds, it could also explain why fluctuations are naturally excited in the time-accurate simulations. Although not 
shown here, similar spikes in growth also exist at other spanwise locations for z < 6 mm. 

The spectral content of the time traces at the probe point PI for the flow under condition III is shown in Fig. 25 for 
comparison. In contrast to condition II, more than one peak is present. Four frequencies, 10, 20, 86, and 128 kHz 
stand out. Some of the other peaks appear to be harmonics of these outstanding modes. It is also interesting to note 
that more higher-frequency modes are present in this case. Figure 26 depicts the temperature rms amplitudes for all 
grid points on the symmetry plane. Compared with Fig. 24 for condition II, there are three distinct peaks at x = -4.3, 
-2.6, -2.2 mm in the front-side separation region. Similar decay followed by growth trend can also be observed for 
this case. Temperature rms amplitudes at four streamwise locations are plotted in Fig. 27 for the flow under 
condition III. The disturbance structure at x = -2 mm appears to be more complex than that under condition II, as 
shown in Fig. 23, perhaps due to the presence of multiple modes. Qualitatively similar evolution of the rms 
amplitudes downstream as in condition II appears to confirm that similar instability mechanisms may also be present 
for this case. 

The sharp growth discussed above was not observed in the flow under condition I or in other calculations for a 
slender cylindrical roughness element configuration investigated in the Purdue University quiet tunnel [15]. Our 
previous 2D investigations also suggest that no spontaneous instability exists for 2D, rectangular roughness element 
under supersonic free-stream conditions [10]. Thus the possible absolute instability observed under conditions II 
and III, if confirmed, can only take place for 3D roughness elements. The additional dimension in the flow allows 
doubly inflectional profiles to be present. This highly unstable profile along with the presence of spanwise modes in 
3D calculations appears to facilitate the spontaneous instability observed here. Since condition I has a smaller 
roughness to boundary-layer thickness ratio (less than 1) and smaller Re k , it is also likely that some threshold values 
in these parameters must be satisfied in order to have absolute instability. Compared to the spontaneous instability 
occurring in a low-speed, 2D rectangular roughness element for which strong instability results in vortex shedding 
right at the roughness site, the current 3D cases have a rapidly decaying region right after the generation of large 
disturbances (see Figs. 24 and 26). This rapidly stabilizing effect may be attributed to the strong favorable pressure 
gradient associated with the supersonic expansion taking place when the flow passes over the edge of the top 
cylinder surface or beyond the center of the roughness element from the cylinder side face. Such stabilizing 
influence apparently offsets the effect of the highly perturbed mode and as a result, no immediate vortex shedding is 
evident around the cylinder. However, the disturbances can undergo a second bout of growth associated with the 
convective instability of the streak and horse-shoe vortices. The possible absolute instability mechanism upstream 
of the roughness plays the role of frequency selector. For the flow under condition II, the dominant 40 kHz mode is 
an outcome of this mechanism. This mode, once excited, persists throughout the whole flowfield and remains 
dominant in the wake region compared to other higher-frequency modes predicted by the secondary instability 
analysis. However, as far as transition is concerned, the amplifying process in the wake may determine how fast 
these disturbances are amplified and whether transition would occur in the wake or not. Condition III appears to 
lead to highly perturbed wake flow. Higher grid resolution computations are under way to determine whether 
transitional flow is present sufficiently far downstream. 

In Danehy’s wind tunnel experiments, it was observed that strong disturbances in the form of corkscrew patterns and 
large disturbances are present along both sides of the cylinder for the flow under condition III (see Fig. 13(d)) [11]. 
The PLIF seeding process favors the visualization of the horse-shoe vortices since it is more difficult for the seeding 
to enter the mushroom streak in the wake. These unstable flow structures in the wake also suggest possible 
spanwise side-swing type of instability for the horse-shoe vortices, as predicted in the present investigation. At this 
2 mm cylinder height, no processed data for conditions I and II is yet available and, therefore, the transition 
mechanisms at these conditions remain unknown at present. 
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Figure 8. Numerical Schlieren at symmetry and Figure 9. Numerical Schlieren pictures at y = 1.5 mm 

exit planes: (a) condition I (b) condition II plane for (a) condition I (b) condition II (c) condition 

(c) condition III III 

Table 3. Coordinates and locations of the probe points for the DTO trip 


Probe 

X 

y 

z 

Location 

D1 

-0.002 

0.0012 

-0.002 

Front side separation region 

D2 

0.002 

0.0006 

0.002 

Right behind the roughness 

D3 

0.005 

0.0006 

0.003 

Near wake 

D4 

0.015 

0.0006 

0.0 

Middle wake 

D5 

0.047 

0.0004 

-0.002 

Far wake 


12 

American Institute of Aeronautics and Astronautics 






Figure 10. Surface limiting streamline patterns and 
vorticity magnitude contours at the exit plane for: (a) 
condition I (b) condition II (c) condition III 


Figure 11. Normalized surface temperature gradient 
distribution for: (a) condition I (b) condition II (c) 
condition III 
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Figure 12. Isosurfaces of vorticity magnitude color 
coded by temperature for: (a) condition I (b) condition 
II (c) condition III 



(d) 

Figure 13. Vorticity magnitude contours at the y = 0.5 
mm plane for: (a) condition I (b) condition II (c) 
condition III (d) PLIF image at a similar height under 
condition III from ref. [11] 
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Figure 14. Vorticity magnitude contours at four 
downstream locations for flow under condition I 
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Figure 16. Vorticity magnitude contours at four 
downstream locations for flow under condition III 
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Figure 15. Vorticity magnitude at four downstream 
locations for flow under condition II 



z 


Figure 17. Streamwise velocity contours at x = 40 
mm under condition II, showing wake streak and 
the horse-shoe vortex 
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Figure 18. Symmetric unstable modes for the streak at x = 40 mm under condition II obtained by 2D 
eigenvalue analysis showing (a) growth rates of various unstable modes and (b) temperature eigenfunctions 
for modes 1 and 2 



Frequency (Hz) Z (mm) 

(a) (b) 

Figure 19. Anti-symmetric unstable modes for the horse-shoe vortex at x = 40 mm under condition II 
obtained by 2D eigenvalue analysis showing (a) growth rates of the most unstable modes and (b) temperature 
eigenfunctions for two frequencies 
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(a) (b) 

Figure 20. Streamwise velocity and temperature time traces at probe P6 with a coordinate of (41.6, 1.9, 5.9) 
and Fourier spectrum content for temperature perturbations at P5 and P6 under condition II: (a) time traces 
(b) Fourier spectrum 



(a) (b) 


Figure 21. Streamwise velocity and temperature rms amplitude contours at x = 40 mm computed by time- 
accurate simulations with self-excited disturbances for condition II (a) streamwise velocity (b) temperature 



Figure 22. Spectral content at probe point PI 
under condition II 


■ i i i ■ 


u: -0.2 0.0 0.3 0 6 0.8 



(a) 

■ MB 


urms: 0.00 0.03 0.06 0.09 



A 
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Figure 23. Streamwise velocity and rms 
amplitude contours at x = -2 mm under 
condition II (a) u velocity (b) rms u amplitude 
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Figure 24. Temperature rms distributions along 
the streamwise direction for all grid points at 
the symmetry plane (z = 0) under condition II 
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Figure 25. Spectral content at probe point PI 
under condition III 
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Figure 26. Temperature rms distributions along 
the streamwise direction for all grid points at 
the symmetry plane (z = 0) under condition III 



0 5 10 15 20 25 

Z 

Figure 27. Contours of rms amplitude of 
temperature fluctuations at four streamwise 
locations under condition III 


The Orbiter DTP Trip 

A scaled-down DTO trip model was tested in NASA Langley’s Mach 10 tunnel with the PLIF technique recently 
[12]. The trip height is about 1mm. Therefore the value of k/S is about 0.65 and 1.25 under conditions II and III, 
respectively. The roughness Reynolds number Re k is estimated to be around 1,300 and 4,630, respectively. The k/5 
value for the DTO trip under condition III is about the same as the cylindrical roughness element under condition II 
but with a smaller Re k value. It is expected that the DTO trip under condition III should have similar yet slightly 
weaker instability characteristics than the cylinder case under condition II. The flow under condition II for this trip 
is relatively more stable and should remain steady. This conjecture is confirmed in the PLIF experiments that 
indeed the flow under condition II is laminar and stable while that under condition III leads to large disturbances 


18 

American Institute of Aeronautics and Astronautics 



downstream [12]. The trip configuration and a tetrahedral mesh are shown in Fig. 28. There are about 43.6 million 
tetrahedral elements in this mesh. The grid distribution is still sub -optimal, mainly due to the over-clustering near 
the underlying block boundaries of the structured domains. Grid improvement and more comprehensive studies are 
still ongoing. Preliminary solutions for the DTO trip under wind tunnel condition II are discussed in this paper. 
Solutions for condition III will be presented after the grid refinement study is completed. 

Five probe locations listed in Table 3 were used to monitor the solution evolution during the time-accurate 
calculations. The resulting time histories of streamwise velocity and temperature are shown in Fig. 29. The 
solutions appear to converge to a steady state and no instability wave was observed. Since the trip height to the 
boundary layer thickness ratio and the roughness Reynolds number are both small, a converged stable solution is to 
be expected. It remains to be seen whether a more refined and better distributed grid would give rise to the self- 
excited wake instability for condition III. Figures 30 -32 show the converged solution by plotting numerical 
Schlieren, surface limiting streamline pattern, and the temperature and vorticity magnitude contours on four cross 
planes at four streamwise locations. The presence of the trip appears to have a stronger effect on the leeward side 
than the windward side, as expected. The flow discontinuity coming off the trip also appears to interact with the 
free-stream shock. A streak of higher temperature is also evident in the trip wake. Even though no self-excited 
disturbances were observed, convective instability is still possible. 



(d) 


Figure 28. The scaled DTO trip wind tunnel model mounted on a flap plate investigated at Langley 
Research Center and the 43.6 million element tetrahedral mesh : (a) model (b) overall mesh (c) center 
and exit planes (d) close-up view near the trip 
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Figure 29. Time histories of streamwise velocity and temperature for time-accurate computations of 
the DTO trip under condition II: (a) u velocity (b) temperature 



Figure 30. Numerical Schlieren picture at the 
center and exit plane under condition II, showing 
shock structures around the DTO trip 



Figure 31. Computed surface limiting streamline 
pattern for the DTO trip under condition II 
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( a ) (b) 

Figure 32. Converged temperature and vorticity magnitude contours at four streamwise locations (in 
mm) for the DTO trip under condition II: (a) temperature (b) vorticity 
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IV. Summary 


The unstructured-mesh, space-time CESE method is used to simulate time-accurate development of the wake flow 
behind an isolated roughness element in a hypersonic boundary layer. Two roughness shapes in the form of a 
cylindrical trip and a scaled version of the DTO trip used in a recent flight experiment have been investigated. The 
tetrahedral mesh with 33.6 million elements has been shown to adequately resolve both the mushroom-like 
centerline streak as well as the horse-shoe vortices wrapping around the cylindrical roughness element. Three wind 
tunnel conditions pertinent to the NASA Langley Mach 10 wind tunnel experiments were computed. The post-shock 
Mach numbers are 6.15 and 4.16. The roughness height to boundary layer thickness ratio varies from 0.8, 1.3, to 
2.5, corresponding to a roughness Reynolds number Re k of about 6,000 - 13,000. It was found that at the smallest 
height ratio, the flow approaches a steady state with a well-formed streak and horse-shoe vortices downstream of the 
roughness. Nevertheless, strong flow unsteadiness in the form of instability waves is present for the two larger 
values of k/S. In particular, for the intermediate value of k/S, instability waves with a dominant frequency of 40 kHz 
were found. This flow instability causes the streak and horse-shoe vortices to oscillate downstream. Spectral 
analyses and the computed evolution of disturbance rms amplitude suggest that a possible absolute instability may 
take place in the front-side separation region where doubly inflectional mean profiles exist on the cross -stream 
plane. Instability waves originating from this region propagate downstream and cause oscillatory motions of the 
streak and the horse-shoe vortices. A similar instability mechanism but with a broader Fourier spectrum also takes 
place for the largest value of k/S. Stronger instability waves both upstream and downstream of the roughness at this 
relatively larger roughness height eventually lead to early vortex breakdown. Detailed flow visualization of these 
processes is documented. Preliminary results are also presented for the wind tunnel model of the DTO trip with a 
height to boundary-layer thickness ratio of about 0.65 and Re k of about 1,300. Time-accurate simulations eventually 
lead to a steady state for the DTI trip. 
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